## Pelc 2011 IO

library(foreign) 
library(Amelia)

## Load original dataset
p2011io <- read.dta("P2011 IO Rep Data.dta")
head(p2011io)
dim(p2011io)

## Drop ID and other index vars
p2011io$reportername <- p2011io$countryname <- p2011io$countrycode <- p2011io$ID <- p2011io$product <- p2011io$tariffyr <- p2011io$membershipyr <- p2011io$appliedyr <- p2011io$implemyrs <- p2011io$joinyr <-  p2011io$controlsyr <- NULL

## Drop derived dummy vars
p2011io$depthdum1 <- p2011io$depthdum2 <-p2011io$post2002 <- p2011io$implemperiod <- p2011io$post2002 <- NULL

## Drop constituent items
p2011io$AHSpre <- p2011io$BNDpost <- p2011io$MFNpost <- p2011io$AHSpost <- p2011io$MFNpreW1 <- p2011io$MFNpostW1 <- p2011io$BNDpreWI <- p2011io$BNDpostWI <- p2011io$AHSpreW1 <- p2011io$AHSpostW1 <- p2011io$AHSpostW <- p2011io$MFNpostW <-p2011io$BNDpostW <-p2011io$AHSpreW <-p2011io$MFNpreW <-p2011io$AHSpostW1 <-p2011io$MFNpostW1 <-p2011io$BNDpostW1 <-p2011io$AHSpreW1 <-p2011io$MFNpreW1 <- p2011io$AHSw <-p2011io$BNDw <-p2011io$MFNw <- NULL

## Drop duplicates
p2011io$implementationyear <- NULL

## Drop vars with no variation
p2011io$defactomember <- NULL

## Drop calculated parameters
p2011io$b1 <- p2011io$b2 <-p2011io$b3 <-p2011io$b4 <-p2011io$b5 <-p2011io$b6 <-p2011io$b7 <-p2011io$b8 <-p2011io$b9 <-p2011io$b10 <- p2011io$b11 <- p2011io$bndfuture <- p2011io$b11 <- p2011io$b11 <-p2011io$b11 <-p2011io$b11 <- p2011io$s_chinad <-p2011io$r_chinad <- p2011io$e_chinad <-p2011io$s_logimp <-p2011io$r_logimp <-p2011io$e_logimp <-p2011io$s_logexp <-p2011io$r_logexp <-p2011io$e_logexp <-p2011io$GAM_mu <-p2011io$s_logful <-p2011io$r_logful <-p2011io$e_logful <-p2011io$s_access <-p2011io$r_access <-p2011io$e_access <-p2011io$s_logGDP <-p2011io$r_logGDP <-p2011io$e_logGDP <-p2011io$s_loggdp <-p2011io$r_loggdp <- p2011io$e_loggdp <-p2011io$s_MFNpre <-p2011io$r_MFNpre <-p2011io$e_MFNpre <-p2011io$s_polity <-p2011io$r_polity <-p2011io$e_polity <-p2011io$s_howLAT <- p2011io$r_howLAT <- p2011io$e_howLAT <- p2011io$CAperc <- p2011io$exportspercent <- p2011io$e_howLAT <- p2011io$depth2percent <- p2011io$depth3percent <- NULL

which( colnames(p2011io)=="_est_t5_1" )
which( colnames(p2011io)=="_est_t5_2" )
p2011io <- p2011io[, -c(81, 82)]

## How many variables? 80: reduction necessary due to large sample size
dim(p2011io)

## Which variables are in analysis and have missing data?
analysis <- as.data.frame(cbind(p2011io$depth2, p2011io$logGDPconst, p2011io$logfullimports, p2011io$logfullexports, p2011io$accessionperiod, p2011io$loggdpcap, p2011io$MFNpre, p2011io$polity3, p2011io$chinadum, p2011io$howLATE, p2011io$reporter, p2011io$depth3, p2011io$Xexp, p2011io$Ximp, p2011io$overhangfinal))

## Create missingness indicator variables
missing <- as.data.frame(cbind(as.integer(complete.cases(p2011io$depth2)), as.integer(complete.cases(p2011io$logGDPconst)), as.integer(complete.cases(p2011io$logfullimports)), as.integer(complete.cases(p2011io$logfullexports)), as.integer(complete.cases(p2011io$accessionperiod)), as.integer(complete.cases(p2011io$loggdpcap)), as.integer(complete.cases(p2011io$MFNpre)), as.integer(complete.cases(p2011io$chinadum)), as.integer(complete.cases(p2011io$polity3)), as.integer(complete.cases(p2011io$howLATE)), as.integer(complete.cases(p2011io$reporter)), as.integer(complete.cases(p2011io$depth3)), as.integer(complete.cases(p2011io$Xexp)), as.integer(complete.cases(p2011io$Ximp)), as.integer(complete.cases(p2011io$overhangfinal))))

dim(analysis)
dim(missing)
apply(missing, 2, sd)

## Remove analysis variables
## p2011io$depth2 <- p2011io$logGDPconst <- p2011io$logfullimports <- p2011io$logfullexports <- p2011io$accessionperiod <- p2011io$loggdpcap <- p2011io$MFNpre <-p2011io$polity3 <- p2011io$chinadum <- p2011io$howLATE <- p2011io$reporter <- p2011io$depth3 <- p2011io$Xexp <- p2011io$Ximp <- p2011io$overhangfinal <- NULL
## head(p2011io)
## var(p2011io)

## Check round(correlations and missing values
round(cor(p2011io$importsvalue000, analysis, use = "pairwise.complete.obs"), 2)
round(cor(p2011io$importsvalue000, missing, use = "pairwise.complete.obs"), 2)
sum(is.na(p2011io$importsvalue000))/nrow(p2011io)*100
## N
p2011io$importsvalue000 <- NULL

round(cor(p2011io$AHS, analysis, use = "pairwise.complete.obs"), 2)
round(cor(p2011io$AHS, missing, use = "pairwise.complete.obs"), 2)
sum(is.na(p2011io$AHS))/nrow(p2011io)*100
## N
p2011io$AHS <- NULL

round(cor(p2011io$BND, analysis, use = "pairwise.complete.obs"), 2)
round(cor(p2011io$BND, missing, use = "pairwise.complete.obs"), 2)
sum(is.na(p2011io$BND))/nrow(p2011io)*100
## N
p2011io$BND <- NULL

round(cor(p2011io$MFN, analysis, use = "pairwise.complete.obs"), 2)
round(cor(p2011io$MFN, missing, use = "pairwise.complete.obs"), 2)
sum(is.na(p2011io$MFN))/nrow(p2011io)*100
## N
p2011io$MFN <- NULL

round(cor(p2011io$PRF, analysis, use = "pairwise.complete.obs"), 2)
round(cor(p2011io$PRF, missing, use = "pairwise.complete.obs"), 2)
sum(is.na(p2011io$PRF))/nrow(p2011io)*100
## N
p2011io$PRF <- NULL

round(cor(p2011io$WTOmember, analysis, use = "pairwise.complete.obs"), 2)
round(cor(p2011io$WTOmember, missing, use = "pairwise.complete.obs"), 2)
sum(is.na(p2011io$WTOmember))/nrow(p2011io)*100
## N
p2011io$WTOmember <- NULL

round(cor(p2011io$PRFpre, analysis, use = "pairwise.complete.obs"), 2)
round(cor(p2011io$PRFpre, missing, use = "pairwise.complete.obs"), 2)
sum(is.na(p2011io$PRFpre))/nrow(p2011io)*100
## N
p2011io$PRFpre <- NULL


round(cor(p2011io$logGDP, analysis, use = "pairwise.complete.obs"), 2)
round(cor(p2011io$logGDP, missing, use = "pairwise.complete.obs"), 2)
sum(is.na(p2011io$logGDP))/nrow(p2011io)*100
## Y

round(cor(p2011io$PRFpost, analysis, use = "pairwise.complete.obs"), 2)
round(cor(p2011io$PRFpost, missing, use = "pairwise.complete.obs"), 2)
sum(is.na(p2011io$PRFpost))/nrow(p2011io)*100
## N
p2011io$PRFpost <- NULL

round(cor(p2011io$depth1, analysis, use = "pairwise.complete.obs"), 2)
round(cor(p2011io$depth1, missing, use = "pairwise.complete.obs"), 2)
sum(is.na(p2011io$depth1))/nrow(p2011io)*100
## N
p2011io$depth1 <- NULL

round(cor(p2011io$depth4, analysis, use = "pairwise.complete.obs"), 2)
round(cor(p2011io$depth4, missing, use = "pairwise.complete.obs"), 2)
sum(is.na(p2011io$depth4))/nrow(p2011io)*100
## N
p2011io$depth4 <- NULL

round(cor(p2011io$currentaccount, analysis, use = "pairwise.complete.obs"), 2)
round(cor(p2011io$currentaccount, missing, use = "pairwise.complete.obs"), 2)
sum(is.na(p2011io$currentaccount))/nrow(p2011io)*100
## Y

round(cor(p2011io$importsvalue000, analysis, use = "pairwise.complete.obs"), 2)
round(cor(p2011io$importsvalue000, missing, use = "pairwise.complete.obs"), 2)
sum(is.na(p2011io$importsvalue000))/nrow(p2011io)*100
## N
p2011io$importsvalue000 <- NULL

round(cor(p2011io$fdi, analysis, use = "pairwise.complete.obs"), 2)
round(cor(p2011io$fdi, missing, use = "pairwise.complete.obs"), 2)
sum(is.na(p2011io$fdi))/nrow(p2011io)*100
## Y

round(cor(p2011io$inflation, analysis, use = "pairwise.complete.obs"), 2)
round(cor(p2011io$inflation, missing, use = "pairwise.complete.obs"), 2)
sum(is.na(p2011io$inflation))/nrow(p2011io)*100
## Y

round(cor(p2011io$aidTOTAL, analysis, use = "pairwise.complete.obs"), 2)
round(cor(p2011io$aidTOTAL, missing, use = "pairwise.complete.obs"), 2)
sum(is.na(p2011io$aidTOTAL))/nrow(p2011io)*100
## Y

round(cor(p2011io$aidEC, analysis, use = "pairwise.complete.obs"), 2)
round(cor(p2011io$aidEC, missing, use = "pairwise.complete.obs"), 2)
sum(is.na(p2011io$aidEC))/nrow(p2011io)*100
## Y

round(cor(p2011io$aidUS, analysis, use = "pairwise.complete.obs"), 2)
round(cor(p2011io$aidUS, missing, use = "pairwise.complete.obs"), 2)
sum(is.na(p2011io$aidUS))/nrow(p2011io)*100
## N
p2011io$aidUS <- NULL

round(cor(p2011io$debt, analysis, use = "pairwise.complete.obs"), 2)
round(cor(p2011io$debt, missing, use = "pairwise.complete.obs"), 2)
sum(is.na(p2011io$debt))/nrow(p2011io)*100
## N
p2011io$debt <- NULL

round(cor(p2011io$GDPconst, analysis, use = "pairwise.complete.obs"), 2)
round(cor(p2011io$GDPconst, missing, use = "pairwise.complete.obs"), 2)
sum(is.na(p2011io$GDPconst))/nrow(p2011io)*100
## N
p2011io$GDPconst <- NULL

round(cor(p2011io$GDPcurr, analysis, use = "pairwise.complete.obs"), 2)
round(cor(p2011io$GDPcurr, missing, use = "pairwise.complete.obs"), 2)
sum(is.na(p2011io$GDPcurr))/nrow(p2011io)*100
## Y

round(cor(p2011io$growth, analysis, use = "pairwise.complete.obs"), 2)
round(cor(p2011io$growth, missing, use = "pairwise.complete.obs"), 2)
sum(is.na(p2011io$growth))/nrow(p2011io)*100
## N
p2011io$growth <- NULL

round(cor(p2011io$imports, analysis, use = "pairwise.complete.obs"), 2)
round(cor(p2011io$imports, missing, use = "pairwise.complete.obs"), 2)
sum(is.na(p2011io$imports))/nrow(p2011io)*100
## N
p2011io$imports <- NULL

round(cor(p2011io$exports, analysis, use = "pairwise.complete.obs"), 2)
round(cor(p2011io$exports, missing, use = "pairwise.complete.obs"), 2)
sum(is.na(p2011io$exports))/nrow(p2011io)*100
## N
p2011io$exports <- NULL

round(cor(p2011io$TOT, analysis, use = "pairwise.complete.obs"), 2)
round(cor(p2011io$TOT, missing, use = "pairwise.complete.obs"), 2)
sum(is.na(p2011io$TOT))/nrow(p2011io)*100
## N
p2011io$TOT <- NULL

round(cor(p2011io$tradepercent, analysis, use = "pairwise.complete.obs"), 2)
round(cor(p2011io$tradepercent, missing, use = "pairwise.complete.obs"), 2)
sum(is.na(p2011io$tradepercent))/nrow(p2011io)*100
## N
p2011io$tradepercent <- NULL

round(cor(p2011io$logCA, analysis, use = "pairwise.complete.obs"), 2)
round(cor(p2011io$logCA, missing, use = "pairwise.complete.obs"), 2)
sum(is.na(p2011io$logCA))/nrow(p2011io)*100
## N
p2011io$logCA <- NULL

round(cor(p2011io$logFDI, analysis, use = "pairwise.complete.obs"), 2)
round(cor(p2011io$logFDI, missing, use = "pairwise.complete.obs"), 2)
sum(is.na(p2011io$logFDI))/nrow(p2011io)*100
## N
p2011io$logFDI <- NULL

round(cor(p2011io$logaidtotal, analysis, use = "pairwise.complete.obs"), 2)
round(cor(p2011io$logaidtotal, missing, use = "pairwise.complete.obs"), 2)
sum(is.na(p2011io$logaidtotal))/nrow(p2011io)*100
## Y

round(cor(p2011io$logaidEC, analysis, use = "pairwise.complete.obs"), 2)
round(cor(p2011io$logaidEC, missing, use = "pairwise.complete.obs"), 2)
sum(is.na(p2011io$logaidEC))/nrow(p2011io)*100
## N
p2011io$logaidEC <- NULL

round(cor(p2011io$logaidUS, analysis, use = "pairwise.complete.obs"), 2)
round(cor(p2011io$logaidUS, missing, use = "pairwise.complete.obs"), 2)
sum(is.na(p2011io$logaidUS))/nrow(p2011io)*100
## Y

round(cor(p2011io$logimports, analysis, use = "pairwise.complete.obs"), 2)
round(cor(p2011io$logimports, missing, use = "pairwise.complete.obs"), 2)
sum(is.na(p2011io$logimports))/nrow(p2011io)*100
## N
p2011io$logimports <- NULL

round(cor(p2011io$logexports, analysis, use = "pairwise.complete.obs"), 2)
round(cor(p2011io$logexports, missing, use = "pairwise.complete.obs"), 2)
sum(is.na(p2011io$logexports))/nrow(p2011io)*100
## N
p2011io$logexports <- NULL

round(cor(p2011io$expprop, analysis, use = "pairwise.complete.obs"), 2)
round(cor(p2011io$expprop, missing, use = "pairwise.complete.obs"), 2)
sum(is.na(p2011io$expprop))/nrow(p2011io)*100
## N
p2011io$expprop <- NULL

round(cor(p2011io$impprop, analysis, use = "pairwise.complete.obs"), 2)
round(cor(p2011io$impprop, missing, use = "pairwise.complete.obs"), 2)
sum(is.na(p2011io$impprop))/nrow(p2011io)*100
## Y

round(cor(p2011io$logcurrentaccount, analysis, use = "pairwise.complete.obs"), 2)
round(cor(p2011io$logcurrentaccount, missing, use = "pairwise.complete.obs"), 2)
sum(is.na(p2011io$logcurrentaccount))/nrow(p2011io)*100
## N
p2011io$logcurrentaccount <- NULL

round(cor(p2011io$ptascore, analysis, use = "pairwise.complete.obs"), 2)
round(cor(p2011io$ptascore, missing, use = "pairwise.complete.obs"), 2)
sum(is.na(p2011io$ptascore))/nrow(p2011io)*100
## N
p2011io$ptascore <- NULL

round(cor(p2011io$openness, analysis, use = "pairwise.complete.obs"), 2)
round(cor(p2011io$openness, missing, use = "pairwise.complete.obs"), 2)
sum(is.na(p2011io$openness))/nrow(p2011io)*100
## Y

round(cor(p2011io$boundLINE, analysis, use = "pairwise.complete.obs"), 2)
round(cor(p2011io$boundLINE, missing, use = "pairwise.complete.obs"), 2)
sum(is.na(p2011io$boundLINE))/nrow(p2011io)*100
## N
p2011io$importsvalue000 <- NULL

round(cor(p2011io$overhangPOST, analysis, use = "pairwise.complete.obs"), 2)
round(cor(p2011io$overhangPOST, missing, use = "pairwise.complete.obs"), 2)
sum(is.na(p2011io$overhangPOST))/nrow(p2011io)*100
## Y

round(cor(p2011io$depthW1, analysis, use = "pairwise.complete.obs"), 2)
round(cor(p2011io$depthW1, missing, use = "pairwise.complete.obs"), 2)
sum(is.na(p2011io$depthW1))/nrow(p2011io)*100
## N
p2011io$depthW1 <- NULL

round(cor(p2011io$depthW2, analysis, use = "pairwise.complete.obs"), 2)
round(cor(p2011io$depthW2, missing, use = "pairwise.complete.obs"), 2)
sum(is.na(p2011io$depthW2))/nrow(p2011io)*100
## N
p2011io$depthW2 <- NULL

round(cor(p2011io$depthW3, analysis, use = "pairwise.complete.obs"), 2)
round(cor(p2011io$depthW3, missing, use = "pairwise.complete.obs"), 2)
sum(is.na(p2011io$depthW3))/nrow(p2011io)*100
## Y

round(cor(p2011io$depthW4, analysis, use = "pairwise.complete.obs"), 2)
round(cor(p2011io$depthW4, missing, use = "pairwise.complete.obs"), 2)
sum(is.na(p2011io$depthW4))/nrow(p2011io)*100
## Y

round(cor(p2011io$BNDpostviet, analysis, use = "pairwise.complete.obs"), 2)
round(cor(p2011io$BNDpostviet, missing, use = "pairwise.complete.obs"), 2)
sum(is.na(p2011io$BNDpostviet))/nrow(p2011io)*100
## N
p2011io$BNDpostviet <- NULL

round(cor(p2011io$BNDpostviet1, analysis, use = "pairwise.complete.obs"), 2)
round(cor(p2011io$BNDpostviet1, missing, use = "pairwise.complete.obs"), 2)
sum(is.na(p2011io$BNDpostviet1))/nrow(p2011io)*100
## N
p2011io$BNDpostviet1 <- NULL

round(cor(p2011io$BNDpostukr, analysis, use = "pairwise.complete.obs"), 2)
round(cor(p2011io$BNDpostukr, missing, use = "pairwise.complete.obs"), 2)
sum(is.na(p2011io$BNDpostukr))/nrow(p2011io)*100
## N
p2011io$BNDpostukr <- NULL

round(cor(p2011io$BNDpostukr1, analysis, use = "pairwise.complete.obs"), 2)
round(cor(p2011io$BNDpostukr1, missing, use = "pairwise.complete.obs"), 2)
sum(is.na(p2011io$BNDpostukr1))/nrow(p2011io)*100
## N
p2011io$BNDpostukr1 <- NULL

round(cor(p2011io$logimportsprod, analysis, use = "pairwise.complete.obs"), 2)
round(cor(p2011io$logimportsprod, missing, use = "pairwise.complete.obs"), 2)
sum(is.na(p2011io$logimportsprod))/nrow(p2011io)*100
## N
p2011io$logimportsprod <- NULL

round(cor(p2011io$policyspace, analysis, use = "pairwise.complete.obs"), 2)
round(cor(p2011io$policyspace, missing, use = "pairwise.complete.obs"), 2)
sum(is.na(p2011io$policyspace))/nrow(p2011io)*100
## N
p2011io$policyspace <- NULL

round(cor(p2011io$pop, analysis, use = "pairwise.complete.obs"), 2)
round(cor(p2011io$pop, missing, use = "pairwise.complete.obs"), 2)
sum(is.na(p2011io$pop))/nrow(p2011io)*100
## N
p2011io$pop <- NULL

round(cor(p2011io$gdpcap, analysis, use = "pairwise.complete.obs"), 2)
round(cor(p2011io$gdpcap, missing, use = "pairwise.complete.obs"), 2)
sum(is.na(p2011io$gdpcap))/nrow(p2011io)*100
## N
p2011io$gdpcap <- NULL

round(cor(p2011io$gdppercapitagrowthannual, analysis, use = "pairwise.complete.obs"), 2)
round(cor(p2011io$gdppercapitagrowthannual, missing, use = "pairwise.complete.obs"), 2)
sum(is.na(p2011io$gdppercapitagrowthannual))/nrow(p2011io)*100
## N
p2011io$gdppercapitagrowthannual <- NULL

round(cor(p2011io$gnipppcurrentinternational, analysis, use = "pairwise.complete.obs"), 2)
round(cor(p2011io$gnipppcurrentinternational, missing, use = "pairwise.complete.obs"), 2)
sum(is.na(p2011io$gnipppcurrentinternational))/nrow(p2011io)*100
## Y

round(cor(p2011io$LDC, analysis, use = "pairwise.complete.obs"), 2)
round(cor(p2011io$LDC, missing, use = "pairwise.complete.obs"), 2)
sum(is.na(p2011io$LDC))/nrow(p2011io)*100
## N
p2011io$LDC <- NULL

round(cor(p2011io$loggdpcap2, analysis, use = "pairwise.complete.obs"), 2)
round(cor(p2011io$loggdpcap2, missing, use = "pairwise.complete.obs"), 2)
sum(is.na(p2011io$loggdpcap2))/nrow(p2011io)*100
## Y

round(cor(p2011io$demXimports, analysis, use = "pairwise.complete.obs"), 2)
round(cor(p2011io$demXimports, missing, use = "pairwise.complete.obs"), 2)
sum(is.na(p2011io$demXimports))/nrow(p2011io)*100
## Y

round(cor(p2011io$fullimports, analysis, use = "pairwise.complete.obs"), 2)
round(cor(p2011io$fullimports, missing, use = "pairwise.complete.obs"), 2)
sum(is.na(p2011io$fullimports))/nrow(p2011io)*100
## N
p2011io$importsvalue000 <- NULL

round(cor(p2011io$fullexports, analysis, use = "pairwise.complete.obs"), 2)
round(cor(p2011io$fullexports, missing, use = "pairwise.complete.obs"), 2)
sum(is.na(p2011io$fullexports))/nrow(p2011io)*100
## N
p2011io$fullexports <- NULL

round(cor(p2011io$prefmargin, analysis, use = "pairwise.complete.obs"), 2)
round(cor(p2011io$prefmargin, missing, use = "pairwise.complete.obs"), 2)
sum(is.na(p2011io$prefmargin))/nrow(p2011io)*100
## N
p2011io$prefmargin <- NULL

round(cor(p2011io$xrcomp, analysis, use = "pairwise.complete.obs"), 2)
round(cor(p2011io$xrcomp, missing, use = "pairwise.complete.obs"), 2)
sum(is.na(p2011io$xrcomp))/nrow(p2011io)*100
## Y

round(cor(p2011io$xropen, analysis, use = "pairwise.complete.obs"), 2)
round(cor(p2011io$xropen, missing, use = "pairwise.complete.obs"), 2)
sum(is.na(p2011io$xropen))/nrow(p2011io)*100
## Y

round(cor(p2011io$xconst, analysis, use = "pairwise.complete.obs"), 2)
round(cor(p2011io$xconst, missing, use = "pairwise.complete.obs"), 2)
sum(is.na(p2011io$xconst))/nrow(p2011io)*100
## Y

round(cor(p2011io$parreg, analysis, use = "pairwise.complete.obs"), 2)
round(cor(p2011io$parreg, missing, use = "pairwise.complete.obs"), 2)
sum(is.na(p2011io$parreg))/nrow(p2011io)*100
## Y

round(cor(p2011io$parcomp, analysis, use = "pairwise.complete.obs"), 2)
round(cor(p2011io$parcomp, missing, use = "pairwise.complete.obs"), 2)
sum(is.na(p2011io$parcomp))/nrow(p2011io)*100
## Y

round(cor(p2011io$effect2002, analysis, use = "pairwise.complete.obs"), 2)
round(cor(p2011io$effect2002, missing, use = "pairwise.complete.obs"), 2)
sum(is.na(p2011io$effect2002))/nrow(p2011io)*100
## N
p2011io$effect2002 <- NULL

## Imputation
case <-c(1:nrow(p2011io))
p2011io <- cbind(p2011io, case)
head(p2011io)
dim(p2011io)
var(p2011io)

## What is average percentage of missing data?
NAs <- function(x) {
    as.vector(apply(x, 2, function(x) length(which(is.na(x)))))
    }
NAs(p2011io)
mean(NAs(p2011io)/nrow(p2011io))*100

## Thus: 16 imputations

set.seed(02138)
p2011io.out <- amelia(p2011io, m = 16, cs = "reporter", empri = 0.01*nrow(p2011io))

write.amelia(obj=p2011io.out, file.stem = "P2011 IO Imp Data", format = "dta", separate = FALSE)